Setup

Notes

  • Cell suspensions of meningeal dura enriched by MACS for cells expressing Cd45, Cd31, and Lyve-1
    • Cd45/Ptprc/B220/Ly-5/Lyt-4/T200: transmembrane protein tyrosine phosphatase, located on most haematopoietic cells, positive selection of leukocytes
    • Cd31/Pecam1: adhesion and signaling receptor that is expressed on endothelial and hematopoietic cells, positive selection of endothelial cells
    • Lyve1/1200012G08Rik/Lyve-1/Xlkd1/lymphatic vessel endothelial HA receptor-1: found primarily on lymphatic endothelial cells

Libraries

library(dplyr)          # ungroup()
library(ggtree)         # BuildClusterTree()
library(gridExtra)      # grid.arrange()
library(gtools)         # smartbind()
library(parallel)       # detectCores()
library(plotly)         # plot_ly()
library(Seurat)         # Idents()
library(SeuratWrappers) # RunPrestoAll()
library(ShinyCell)      # createConfig()
library(tidyr)          # %>%

Variables

out <- "../../results/all_clusters/"
sample_order <- c("E3.M","E3.F","E4.M","E4.F")
sample_colors <- c("#26946A","#1814A1","#EAC941","#DF5F00")
sample_order2 <- c("Male E3", "Male E4", "Female E3", "Female E4")
isoform_order <- c("E4","E3")
isoform_colors <- c("darkgray","cornflowerblue")
sex_order <- c("Male","Female")
sex_colors <- c("darkgray","purple")

Read data

mouse.unannotated <- readRDS("../../rObjects/unannotated_obj.rds")
Idents(mouse.unannotated) <- "seurat_clusters"
DefaultAssay(mouse.unannotated) <- "RNA"
mouse.unannotated <- NormalizeData(mouse.unannotated)
## Normalizing layer: counts
mouse.unannotated
## An object of class Seurat 
## 41866 features across 16965 samples within 2 assays 
## Active assay: RNA (21089 features, 0 variable features)
##  2 layers present: data, counts
##  1 other assay present: SCT
##  3 dimensional reductions calculated: pca, harmony, umap

Unannotated QC

UMAP

cluster_colors <- c("red4","red","orange","gold","yellow","lightgreen","chartreuse3","darkgreen",
                    "cyan","steelblue","blue","deeppink","pink","salmon","purple3","orchid","tan","chocolate4",
                    "gray30","gray70","black","aquamarine")
u1 <- DimPlot(object = mouse.unannotated,
              reduction = "umap",
              shuffle = TRUE,
              repel = TRUE,
              raster = FALSE,
              cols = cluster_colors,
              label = TRUE)
u1

u2 <- DimPlot(object = mouse.unannotated,
              reduction = "umap",
              shuffle = TRUE,
              repel = TRUE,
              dims = c(2,3),
              cols = cluster_colors,
              label = TRUE)
u2

Feature plots

# UMAP percent.mt
FeaturePlot(mouse.unannotated,
            reduction = "umap", 
            features = "percent.mt") + 
  scale_colour_gradientn(colours = c("blue","lightblue","yellow","orange","red"))

# UMAP percent.ribo
FeaturePlot(mouse.unannotated,
            reduction = "umap", 
            features = "percent.ribo.protein") + 
  scale_colour_gradientn(colours = c("blue","lightblue","yellow","orange","red"))

# UMAP percent.hb
FeaturePlot(mouse.unannotated,
            reduction = "umap", 
            features = "percent.hb") + 
  scale_colour_gradientn(colours = c("blue","lightblue","yellow","orange","red"))

# UMAP nCount
FeaturePlot(mouse.unannotated,
            reduction = "umap", 
            features = "nCount_RNA") + 
  scale_colour_gradientn(colours = c("blue","lightblue","yellow","orange","red"))

# UMAP nFeature
FeaturePlot(mouse.unannotated,
            reduction = "umap", 
            features = "nFeature_RNA") + 
  scale_colour_gradientn(colours = c("blue","lightblue","yellow","orange","red"))

# UMAP cell.complexity
FeaturePlot(mouse.unannotated,
            reduction = "umap", 
            features = "cell.complexity") + 
  scale_colour_gradientn(colours = c("blue","lightblue","yellow","orange","red"))

Violins

VlnPlot(mouse.unannotated,
        features = "nCount_RNA",
        split.by = "seurat_clusters")
## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##       
## This message will be shown once per session.

VlnPlot(mouse.unannotated,
        features = "nFeature_RNA",
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "cell.complexity",
        split.by = "seurat_clusters")

Cells per cluster

# Cells per sample per cluster
sample_ncells <- FetchData(mouse.unannotated, 
                     vars = c("ident", "sample")) %>%
  dplyr::count(ident,sample) %>%
  tidyr::spread(ident, n)
sample_ncells
##   sample   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16 17
## 1   E3.M 598 529 228 220 267 254 189 282 344 196 118 110  84 105  94  49  29 94
## 2   E3.F 288 578 302 326 199 229 184 285  27 116 178  60  65  31  40  21  40 58
## 3   E4.M 581 381 238 218 232 242 228 304 313 238 165 113  63  66  74  55  37 48
## 4   E4.F 818 769 596 425 366 312 419 109 289 337 295 151 155 142 122 182 171 70
##   18 19 20 21
## 1 41 64 59 18
## 2 70  7 43  8
## 3 62 80 62 21
## 4 82 93 59 55
# Cells per isoform per cluster
isoform_ncells <- FetchData(mouse.unannotated, 
                     vars = c("ident", "isoform")) %>%
  dplyr::count(ident,isoform) %>%
  tidyr::spread(ident, n)
isoform_ncells
##   isoform    0    1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16
## 1      E4 1399 1150 834 643 598 554 647 413 602 575 460 264 218 208 196 237 208
## 2      E3  886 1107 530 546 466 483 373 567 371 312 296 170 149 136 134  70  69
##    17  18  19  20 21
## 1 118 144 173 121 76
## 2 152 111  71 102 26
# Cells per sex per cluster
sex_ncells <- FetchData(mouse.unannotated, 
                     vars = c("ident", "sex")) %>%
  dplyr::count(ident,sex) %>%
  tidyr::spread(ident, n)
sex_ncells
##      sex    0    1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16
## 1   Male 1179  910 466 438 499 496 417 586 657 434 283 223 147 171 168 104  66
## 2 Female 1106 1347 898 751 565 541 603 394 316 453 473 211 220 173 162 203 211
##    17  18  19  20 21
## 1 142 103 144 121 39
## 2 128 152 100 102 63

Gene histogram

# User params
goi <- "Malat1"
threshold <- 1

# Subset data
log2.threshold <- log2(threshold + 0.01)
counts.df <- FetchData(mouse.unannotated, vars = goi)
colnames(counts.df) <- "counts"
log2.counts.df <- log2(counts.df + 0.01)

# Histogram
title <- paste0(goi, "\nnCount_RNA > ", threshold)
hist1 <- ggplot(counts.df, aes(x = counts)) + 
  geom_histogram(bins = 100, fill = "gray", color = "black") + 
  labs(title = title, x=NULL, y=NULL) +
  xlab(paste0(goi, " nCount_RNA")) + ylab("# of Samples") + theme_bw() +
  geom_vline(xintercept = threshold, col = "blue") +
  annotate("rect", 
           xmin = -Inf,
           xmax = threshold, 
           ymin = 0, 
           ymax=Inf, 
           alpha=0.2, 
           fill="chocolate4") +
  annotate("rect", 
           xmin = threshold,
           xmax = Inf, 
           ymin = 0, 
           ymax=Inf, 
           alpha=0.2, 
           fill="deepskyblue")

# Histogram log transformed
hist2 <- ggplot(log2.counts.df, aes(x = counts)) + 
  geom_histogram(bins = 100, fill = "gray", color = "black") + 
  labs(title = title, x=NULL, y=NULL) +
  xlab(paste0("Log2(",goi, " nCount_RNA)")) + ylab("# of Samples") + theme_bw() +
  geom_vline(xintercept = log2.threshold, col = "blue") +
  annotate("rect", 
           xmin = -Inf,
           xmax = log2.threshold, 
           ymin = 0, 
           ymax=Inf, 
           alpha=0.2, 
           fill="chocolate4") +
  annotate("rect", 
           xmin = log2.threshold,
           xmax = Inf, 
           ymin = 0, 
           ymax=Inf, 
           alpha=0.2, 
           fill="deepskyblue")

# plot
plots1 <- list(hist1,hist2)
layout1 <- rbind(c(1),c(2))
grid1 <- grid.arrange(grobs = plots1, layout_matrix = layout1)

Percent gene

# user define variable
goi <- "Malat1"

# Extract counts data
DefaultAssay(mouse.unannotated) <- "RNA"
Idents(mouse.unannotated) <- "seurat_clusters"
geneData <- FetchData(mouse.unannotated,
                      vars = goi,
                      slot = "counts")
## Warning: The `slot` argument of `FetchData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
geneData <- geneData > 1
table(geneData)
## geneData
## FALSE  TRUE 
##  1248 15717
mouse.unannotated$Expression <- geneData

FetchData(mouse.unannotated,
          vars = c("ident", "Expression")) %>%
  dplyr::count(ident, Expression) %>%
  tidyr::spread(ident, n)
##   Expression    0    1    2    3    4    5   6   7   8   9  10  11  12  13  14
## 1      FALSE    2    3    4    2   43    1 690   9  NA  NA  NA  NA 285  13 122
## 2       TRUE 2283 2254 1360 1187 1021 1036 330 971 973 887 756 434  82 331 208
##    15  16  17  18  19  20  21
## 1   2  65   1   4   1   1  NA
## 2 305 212 269 251 243 222 102
# Plot
mouse.unannotated@meta.data %>%
  group_by(seurat_clusters, Expression) %>%
  dplyr::count() %>%
  group_by(seurat_clusters) %>%
  dplyr::mutate(percent = 100*n/sum(n)) %>%
  ungroup() %>%
  ggplot(aes(x=seurat_clusters,y=percent, fill=Expression)) +
  geom_col() +
  ggtitle(paste0("Percentage of cells with > 1 counts for ", goi)) +
  theme(axis.text.x = element_text(angle = 90))

Cluster tree

  • Cluster trees are helpful in deciding what clusters to merge.
mouse.unannotated <- BuildClusterTree(object = mouse.unannotated,
                                     dims = 1:15,
                                     reorder = FALSE,
                                     reorder.numeric = FALSE)
tree <- mouse.unannotated@tools$BuildClusterTree
tree$tip.label <- paste0(tree$tip.label)

ggtree::ggtree(tree, aes(x, y)) +
  scale_y_reverse() +
  ggtree::geom_tree() +
  ggtree::theme_tree() +
  ggtree::geom_tiplab(offset = 1) +
  ggtree::geom_tippoint(color = cluster_colors[1:length(tree$tip.label)], shape = 16, size = 5) +
  coord_cartesian(clip = 'off') +
  theme(plot.margin = unit(c(0,2.5,0,0), 'cm'))

Potential Markers

MACS markers

  • Ptprc - located on most haematopoietic cells, positive selection of leukocytes
  • Pecam1 - adhesion and signaling receptor that is expressed on endothelial and hematopoietic cells, positive selection of endothelial cells
  • Lyve1 - found primarily on lymphatic endothelial cells
VlnPlot(mouse.unannotated,
        features = "Ptprc",
        cols = cluster_colors,
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "Lyve1",
        cols = cluster_colors,
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "Pecam1",
        cols = cluster_colors,
        split.by = "seurat_clusters")

B-cells / Plasma cells

  • Cd19: expressed in B-cells and follicular dendritic cells
  • Fcrla: a B-cell specific protein in mice
  • Cd79a & Cd79b: together form BCR complex
  • Sdc1: Plasma cells
VlnPlot(mouse.unannotated,
        features = "Cd19",
        cols = cluster_colors,
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "Fcrla",
        cols = cluster_colors,
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "Cd79a",
        cols = cluster_colors,
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "Cd79b",
        cols = cluster_colors,
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "Sdc1",
        cols = cluster_colors,
        split.by = "seurat_clusters")

T-cells

  • Trac/Cd3d/Cd3e/Cd3g are components of the TCR
VlnPlot(mouse.unannotated,
        features = "Trac",
        cols = cluster_colors,
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "Cd3e",
        cols = cluster_colors,
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "Cd8a",
        cols = cluster_colors,
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "Cd4",
        cols = cluster_colors,
        split.by = "seurat_clusters")

Endothelial cells

VlnPlot(mouse.unannotated,
        features = "Ly6c1",
        cols = cluster_colors,
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "Flt1",
        cols = cluster_colors,
        split.by = "seurat_clusters")

Fibroblasts

VlnPlot(mouse.unannotated,
        features = "Col1a1",
        cols = cluster_colors,
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "Col1a2",
        cols = cluster_colors,
        split.by = "seurat_clusters")

Mast cells

VlnPlot(mouse.unannotated,
        features = "Fcer1a",
        cols = cluster_colors,
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "Kit", # aka Cd117
        cols = cluster_colors,
        split.by = "seurat_clusters")

Macrophage

VlnPlot(mouse.unannotated,
        features = "C1qa",
        cols = cluster_colors,
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "C1qb",
        cols = cluster_colors,
        split.by = "seurat_clusters")

Monocytes/DCs

  • Cd11c/Itgax: Monocytes & DCs
VlnPlot(mouse.unannotated,
        features = c("Itgax","Cd209a","Flt3","Zbtb46","Ccr2"),
        cols = cluster_colors,
        split.by = "seurat_clusters",
        stack = TRUE,
        flip = TRUE)

Neutrophils

VlnPlot(mouse.unannotated,
        features = "Ly6g",
        cols = cluster_colors,
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "Retnlg",
        cols = cluster_colors,
        split.by = "seurat_clusters")

Pericytes & SMCs

VlnPlot(mouse.unannotated,
        features = "Acta2",
        cols = cluster_colors,
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "Myl9",
        cols = cluster_colors,
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "Rgs5",
        cols = cluster_colors,
        split.by = "seurat_clusters")

Schwann cells

VlnPlot(mouse.unannotated,
        features = "Cdh19",
        cols = cluster_colors,
        split.by = "seurat_clusters")

VlnPlot(mouse.unannotated,
        features = "Mpz",
        cols = cluster_colors,
        split.by = "seurat_clusters")

Sandro’s markers

Idents(mouse.unannotated) <- "seurat_clusters"
goi <- c("Pdgfrb","Acta2","Vwf","Cldn5","Pecam1","Flt4","Prox1","Lyve1","Ptprc",
         "Cd19","Ms4a1","Ighd","Igha","Sdc1","Cd3e","Trbc2","Il7r","Nkg7","Klrb1b","Klrb1c",
         "Gata3","Rora","Itgax","H2-Eb1","Ccr2","Ly6c2","Lyz2","Ly6g","Itgam",
         "Mrc1","Csf1r","Cd38","Mki67","Mcpt4","Ms4a2","Col1a2","Plp1")

v1 <- VlnPlot(mouse.unannotated,
              features = goi,
              cols = cluster_colors,
              split.by = "seurat_clusters",
              flip = TRUE,
              stack = TRUE)
v1

pdf(paste0(out, "markers/sandro_markers_violin_unannotated.pdf"), width = 8, height = 8)
v1
dev.off()
## png 
##   2

Automatically detect markers

# auto find markers
Idents(mouse.unannotated) <- "seurat_clusters"
all.markers <- SeuratWrappers::RunPrestoAll(
  object = mouse.unannotated,
  assay = "RNA",
  slot = "counts",
  only.pos = TRUE
)
## Calculating cluster 0
## Calculating cluster 1
## Calculating cluster 2
## Calculating cluster 3
## Calculating cluster 4
## Calculating cluster 5
## Calculating cluster 6
## Calculating cluster 7
## Calculating cluster 8
## Calculating cluster 9
## Calculating cluster 10
## Calculating cluster 11
## Calculating cluster 12
## Calculating cluster 13
## Calculating cluster 14
## Calculating cluster 15
## Calculating cluster 16
## Calculating cluster 17
## Calculating cluster 18
## Calculating cluster 19
## Calculating cluster 20
## Calculating cluster 21
# filter based on p_val_adj
all.markers <- all.markers[all.markers$p_val_adj < 0.01,]

# subset table to top 2 and top 20 markers
top2 <- Reduce(rbind,
               by(all.markers,
                  all.markers["cluster"],
                  head,
                  n = 2))
top20 <- Reduce(rbind,
               by(all.markers,
                  all.markers["cluster"],
                  head,
                  n = 20))

# plot violin
v1 <- VlnPlot(mouse.unannotated,
        features = top2$gene,
        split.by = "seurat_clusters",
        flip = TRUE,
        stack = TRUE,
        cols = cluster_colors)
v1

# save table
write.table(all.markers, 
            paste0(out, "markers/unannotated_auto_find_markers_adjpval_0.01.tsv"),
            quote = FALSE,
            row.names = FALSE)
# compare 
table(all.markers$cluster)
## 
##    0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
## 2243 3010 3430 2888 1618 3451 1322 1304 2491 1382 1344 4527  570 2207  491  867 
##   16   17   18   19   20   21 
##  536  649 1739 1986 1138  538
# subset based on cluster
cluster0 <- all.markers[all.markers$cluster == 0,]
cluster1 <- all.markers[all.markers$cluster == 1,]
cluster2 <- all.markers[all.markers$cluster == 2,]
cluster3 <- all.markers[all.markers$cluster == 3,]
cluster4 <- all.markers[all.markers$cluster == 4,]
cluster5 <- all.markers[all.markers$cluster == 5,]
cluster6 <- all.markers[all.markers$cluster == 6,]
cluster7 <- all.markers[all.markers$cluster == 7,]
cluster8 <- all.markers[all.markers$cluster == 8,]
cluster9 <- all.markers[all.markers$cluster == 9,]
cluster10 <- all.markers[all.markers$cluster == 10,]
cluster11 <- all.markers[all.markers$cluster == 11,]
cluster12 <- all.markers[all.markers$cluster == 12,]
cluster13 <- all.markers[all.markers$cluster == 13,]
cluster14 <- all.markers[all.markers$cluster == 14,]
cluster15 <- all.markers[all.markers$cluster == 15,]
cluster16 <- all.markers[all.markers$cluster == 16,]
cluster17 <- all.markers[all.markers$cluster == 17,]
cluster18 <- all.markers[all.markers$cluster == 18,]
cluster19 <- all.markers[all.markers$cluster == 19,]
cluster20 <- all.markers[all.markers$cluster == 20,]
cluster21 <- all.markers[all.markers$cluster == 21,]

Cluster Annotations

Cluster 0: BECs

VlnPlot(mouse.unannotated,
        features = cluster0$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 1: Macrophages

VlnPlot(mouse.unannotated,
        features = cluster1$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 2: Fibroblasts

VlnPlot(mouse.unannotated,
        features = cluster2$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 3: Macrophages

VlnPlot(mouse.unannotated,
        features = cluster3$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 4: BECs

VlnPlot(mouse.unannotated,
        features = cluster4$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 5: Dendritic cells

VlnPlot(mouse.unannotated,
        features = cluster5$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 6: Macrophages

VlnPlot(mouse.unannotated,
        features = cluster6$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 7: Neutrophils

VlnPlot(mouse.unannotated,
        features = cluster7$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 8: B cells

VlnPlot(mouse.unannotated,
        features = cluster8$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 9: T and NK cells

VlnPlot(mouse.unannotated,
        features = cluster9$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 10: ILCs

VlnPlot(mouse.unannotated,
        features = cluster10$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 11: T and NK cells

VlnPlot(mouse.unannotated,
        features = cluster11$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 12: Fibroblasts

VlnPlot(mouse.unannotated,
        features = cluster12$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 13: Pericytes and SMCs

VlnPlot(mouse.unannotated,
        features = cluster13$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 14: BECs

VlnPlot(mouse.unannotated,
        features = cluster14$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 15: Schwann cells

VlnPlot(mouse.unannotated,
        features = cluster15$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 16: Fibroblasts

VlnPlot(mouse.unannotated,
        features = cluster16$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 17: Mast cells

VlnPlot(mouse.unannotated,
        features = cluster17$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 18: Myeloid precursors

VlnPlot(mouse.unannotated,
        features = cluster18$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 19: B cells

VlnPlot(mouse.unannotated,
        features = cluster19$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 20: Plasma cells

VlnPlot(mouse.unannotated,
        features = cluster20$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Cluster 21: LECs

VlnPlot(mouse.unannotated,
        features = cluster21$gene[1:20],
        cols = cluster_colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

Assign Identities

# Rename all identities
mouse.annotated <- RenameIdents(object = mouse.unannotated,
                                "0" = "BECs",
                                "1" = "Macrophages",
                                "2" = "Fibroblasts",
                                "3" = "Macrophages",
                                "4" = "BECs",
                                "5" = "Dendritic cells",
                                "6" = "Macrophages",
                                "7" = "Neutrophils",
                                "8" = "B cells",
                                "9" = "T and NK cells",
                                "10" = "ILCs",
                                "11" = "T and NK cells",
                                "12" = "Fibroblasts",
                                "13" = "Pericytes and SMCs",
                                "14" = "BECs",
                                "15" = "Schwann cells",
                                "16" = "Fibroblasts",
                                "17" = "Mast cells",
                                "18" = "Myeloid precursors",
                                "19" = "B cells",
                                "20" = "Plasma cells",
                                "21" = "LECs")

# save idents
mouse.annotated$annotated_clusters <- Idents(mouse.annotated)

# set levels
mouse.annotated$annotated_clusters <- factor(mouse.annotated$annotated_clusters,
                                             levels = c("Pericytes and SMCs",
                                                        "BECs",
                                                        "LECs",
                                                        "B cells",
                                                        "Plasma cells",
                                                        "T and NK cells",
                                                        "ILCs",
                                                        "Dendritic cells",
                                                        "Neutrophils",
                                                        "Macrophages",
                                                        "Myeloid precursors",
                                                        "Mast cells",
                                                        "Fibroblasts",
                                                        "Schwann cells"))

# set ident
Idents(mouse.annotated) <- "annotated_clusters"

Annotated UMAP

# set params
DefaultAssay(mouse.annotated) <- "RNA"
mouse.annotated <- NormalizeData(mouse.annotated)
## Normalizing layer: counts
cluster_colors <- c("#B5B9BA", # Pericytes and SMCs
                    "#3385BB", # BECs
                    "#40BBFF", # LECs
                    "#a5d5a9", # B cells
                    "#5dbd64", # Plasma cells
                    "#1c7e24", # T and NK cells
                    "#F57C7C", # ILCs
                    "#E42622", # Dendritic cells
                    "#FBB268", # Neutrophils
                    "#FE8D19", # Macrophages
                    "#DE9E83", # Myeloid precursors
                    "#A6CEE3", # Mast cells
                    "#9D7BBA", # Fibroblasts
                    "#977899") # Schwann cells

# umap
umap1 <- DimPlot(object = mouse.annotated, 
                reduction = "umap",
                repel = TRUE,
                group.by = "annotated_clusters",
                cols = cluster_colors)
umap1

umap2 <- DimPlot(object = mouse.annotated, 
                reduction = "umap",
                repel = TRUE,
                dims = c(2,3),
                group.by = "annotated_clusters",
                cols = cluster_colors)
umap2

Automatically detect markers

# auto find markers
Idents(mouse.annotated) <- "annotated_clusters"
all.markers <- SeuratWrappers::RunPrestoAll(
  object = mouse.annotated,
  assay = "RNA",
  slot = "counts",
  only.pos = TRUE
)
## Calculating cluster Pericytes and SMCs
## Calculating cluster BECs
## Calculating cluster LECs
## Calculating cluster B cells
## Calculating cluster Plasma cells
## Calculating cluster T and NK cells
## Calculating cluster ILCs
## Calculating cluster Dendritic cells
## Calculating cluster Neutrophils
## Calculating cluster Macrophages
## Calculating cluster Myeloid precursors
## Calculating cluster Mast cells
## Calculating cluster Fibroblasts
## Calculating cluster Schwann cells
# filter based on p_val_adj
all.markers <- all.markers[all.markers$p_val_adj < 0.01,]

# subset table to top 2 and top 20
top2 <- Reduce(rbind,
               by(all.markers,
                  all.markers["cluster"],
                  head,
                  n = 2))
top20 <- Reduce(rbind,
               by(all.markers,
                  all.markers["cluster"],
                  head,
                  n = 20))

# save
write.table(all.markers, 
            paste0(out, "markers/annotated_auto_find_markers_adjpval_0.01.tsv"),
            quote = FALSE,
            row.names = FALSE)

Violin markers

v1 <- VlnPlot(mouse.annotated,
              features = top2$gene,
              split.by = "annotated_clusters",
              flip = TRUE,
              stack = TRUE,
              cols = cluster_colors)
v1

goi <- c("Pdgfrb","Acta2","Vwf","Cldn5","Pecam1","Flt4","Prox1","Lyve1","Ptprc",
         "Cd19","Ms4a1","Ighd","Igha","Sdc1","Cd3e","Trbc2","Il7r","Nkg7","Klrb1b","Klrb1c",
         "Gata3","Rora","Itgax","H2-Eb1","Ccr2","Ly6c2","Lyz2","Ly6g","Itgam",
         "Mrc1","Csf1r","Cd38","Mki67","Mcpt4","Ms4a2","Col1a2","Plp1")
v2 <- VlnPlot(mouse.annotated,
              group.by = "annotated_clusters",
              split.by = "annotated_clusters",
              stack = TRUE,
              flip = TRUE,
              cols = cluster_colors,
              features = goi)
v2

Pseudobulk PCA

  • Each sample will have data points equal to the number of annotated_clusters
# Step 1: Pseudo-bulk the counts based on sample and cell type
pseudo <- AggregateExpression(
    object = mouse.annotated, 
    assays = "RNA", 
    features = rownames(mouse.annotated),
    return.seurat = TRUE, 
    group.by = c("sample", "annotated_clusters")
)
## Centering and scaling data matrix
# Step 2: Normalize the data
pseudo <- NormalizeData(pseudo,
                        normalization.method = "LogNormalize",
                        scale.factor = 10000)
## Normalizing layer: counts
# Step 3: Find variable features
pseudo <- FindVariableFeatures(pseudo,
                               selection.method = "vst", 
                               nfeatures = 2000)
## Finding variable features for layer counts
# Step 4: Scale the data
pseudo <- ScaleData(pseudo, 
                    features = rownames(pseudo))
## Centering and scaling data matrix
# Step 5: Run PCA
pseudo <- RunPCA(pseudo, 
                 features = VariableFeatures(pseudo),
                 npcs = 10)
## PC_ 1 
## Positive:  Dclre1c, Bcl2a1b, Ms4a6b, Fcgr2b, Rgs10, Gpr65, Cd86, Ctsc, Cfp, Gm26542 
##     Clcn5, Tmem8, Zeb2os, Fcgr3, Gpr160, Slc37a2, Nlrc4, Apobec1, Csf1r, P2ry12 
##     Cd300c2, Ccr5, Adgre1, Gm21188, Ifi213, Rbpj, Clec10a, Cd68, Cndp2, C1qb 
## Negative:  Grb10, Plpp3, Fbxl7, Cttn, Ppic, Ghr, Uaca, Smad6, Synpo, Cachd1 
##     Adarb1, Yap1, Serpinh1, Tmtc1, Bgn, Sulf1, Col5a2, Tcaf1, Wwtr1, Pdgfd 
##     Fam20a, Prdm5, Prickle1, Nfib, Fat4, Klf12, Mn1, Nr1d1, Id1, Nid1 
## PC_ 2 
## Positive:  Serpinf1, Slc35g1, Olfml3, Clmp, Frmpd4, Shisa3, Car13, H19, Egfr, Col12a1 
##     Ecrg4, Cdon, Sfrp4, Mfap4, Shank1, Clec11a, Lmx1b, Igdcc4, Cpxm1, Prrx2 
##     Six2, Mdk, Gm29478, Kcns1, Foxd1, C130021I20Rik, Cacna2d2, Slc4a10, Fzd7, Erc2 
## Negative:  Dll4, Rasip1, Exoc3l2, Cd300lg, Arhgef15, Icam2, Mmrn2, Tie1, Cmtm8, Unc45b 
##     Plvap, Tek, Tead4, Ushbp1, Egfl7, Clic5, Sox18, Kdr, Rassf9, Pcdh12 
##     Dipk2b, Arl15, Nova2, Myct1, Gpihbp1, Spns2, Galnt15, Pecam1, Kank3, Btnl9 
## PC_ 3 
## Positive:  Arhgap24, Chchd10, Scd1, Gfra1, Srpk3, Cpm, Chst3, Nrxn1, Fcrla, Sbk1 
##     Gm42836, Spib, Cecr2, 2010309G21Rik, Ctnna3, Gm43388, Ly6d, Cd79a, Bfsp2, Tmem108 
##     Siglecg, Agbl1, Vpreb3, Il9r, B3gnt5, Tnfrsf13c, Ms4a1, H2-Eb2, Pou2af1, Klhl14 
## Negative:  Dab2, Gm4951, Rab3il1, Wwp1, Slco2b1, Ophn1, Stab1, Rnase4, Cmklr1, Ang 
##     Trpv4, F830016B08Rik, Hpgd, Prune2, Serpinb8, Vcam1, Abca9, Gm15635, Mrc1, Ctsb 
##     Cd36, C3ar1, Reps2, 2610203C22Rik, Fcrls, Lgmn, Hfe, Nfxl1, Ms4a7, Msr1 
## PC_ 4 
## Positive:  Dennd3, Erg, Aff3, 2700054A10Rik, Nxn, Gm16046, Cacna1e, Rgs7, Chrm3, Eya1 
##     Kcnj15, Gm16070, Ifitm10, Pik3c2b, Tnfsf10, Fcrla, Ly6d, Cd79a, Dchs1, Myzap 
##     Spib, Siglecg, Trpm6, 2010309G21Rik, Bach2, Pde2a, Gm43388, Pecam1, Gm26760, Dennd5b 
## Negative:  Kcnj12, Itga7, Ctnna3, Axl, Il34, Sncg, Tagln, Mertk, Nrip2, Rasl12 
##     Myom1, Synpo2, Akap6, Cap2, Myh11, Myl9, Dgkb, Adap2, Pln, Gm34829 
##     Jph2, Susd5, Rgs6, Tbx3os1, Casq2, Ppp1r14a, Tpm2, Trarg1, Rgs4, Acta2 
## PC_ 5 
## Positive:  Mmp8, Lcn2, Cd177, Anxa1, Trem3, Dhrs9, B430306N03Rik, 4930438A08Rik, Scrg1, B230208H11Rik 
##     Cebpe, Acod1, Il1rn, Wfdc21, Camp, Slfn4, Fpr2, F5, Arg2, Fpr1 
##     Ngp, Clec4d, A530064D06Rik, Trem1, Padi4, Ly6g, Itgb2l, 9830107B12Rik, Cxcr2, Ltf 
## Negative:  Ptprcap, Gimap7, Ikzf3, Cd2, Gimap3, Fam169b, Gpr174, Bach2, Unc5cl, Ablim1 
##     Aff3, Slamf6, Sidt1, Myo3b, Tox, Rnf43, Lck, Ctla4, A430093F15Rik, Blk 
##     Clnk, Btla, Trbc2, Gm39323, Il2rb, Hdac7, Themis, Gprin3, Ctsw, Zap70
# Step 6: Visualize PCA
pca_colors <- c("firebrick1","cyan","gold","blue","black","forestgreen",
                    "darkorchid1","green","gray","deeppink","chocolate4",
                    "steelblue","pink","orange")
pca <- DimPlot(pseudo, 
        reduction = "pca",
        group.by = "annotated_clusters",
        cols = pca_colors,
        pt.size = 3)
pca

pdf(paste0(out, "clustering_QC/pseudobulk_pca.pdf"), width = 6, height = 4)
pca
dev.off()
## png 
##   2

Shiny App

# create new object
shiny.obj <- mouse.annotated
VariableFeatures(shiny.obj) <- shiny.obj@assays$SCT@var.features

# set default params
DefaultAssay(shiny.obj) <- "RNA"
Idents(shiny.obj) <- "annotated_clusters"

# create config
names <- colnames(shiny.obj@meta.data)
names <- names[c(22,23,2:21)]
sc.config <- createConfig(obj = shiny.obj,
                          meta.to.include = names)

# change wd
setwd(out)

# output shiny app folder
makeShinyApp(obj = shiny.obj, 
             scConf = sc.config, 
             gene.mapping = TRUE,
             shiny.title = "All Clusters")

# manual config edits
sc1conf <- readRDS("shinyApp/sc1conf.rds")
sc1conf[2,4] <- "#B5B9BA|#3385BB|#40BBFF|#a5d5a9|#5dbd64|#1c7e24|#F57C7C|#E42622|#FBB268|#FE8D19|#DE9E83|#A6CEE3|#9D7BBA|#977899"
saveRDS(sc1conf, "shinyApp/sc1conf.rds")